using System;
using System.Reflection;
using System.IO;
using NBody;
using NBody.Analysis;
using OptionParser;
using DataFile;

[assembly: AssemblyTitle ("NBodyNumberDensity")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Calculate the number density N(E, j).")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class NumberDensityOptions : InputOutputGetOptions
{
    [Option ("-b", "Number of bins")]
    public int NumBins;
     
    [Option ("-c", "Adjust for centre of mass")]
    public bool CentreOfMass;
    
    [Option ("-e", "Softening length")]
    public double Softening;
    
    [Option ("-f", "File for potential data")]
    public string PotFile;
    
    [Option ("-B", "Output as BASIC binary image file")]
    public bool Binary;
    
    [Option ("-G", "Output as text, but for Gnuplot")]
    public bool Gnuplot;
    
    [Option ("-M", "Output as text, but for Maple")]
    public bool Maple;
    
    public NumberDensityOptions()
    {
        NumBins = 100;
        Softening = 0.0;
        PotFile = "";
    }
}

public class NBodyNumberDensity
{
    public static void Main(string[] args)
    {
        NumberDensityOptions opts = new NumberDensityOptions();
        opts.ProcessArgs(args, ".numden");
                        
        NBodySystem s = NBodySystem.Read(opts.InFile);
        
        if (opts.CentreOfMass)
            s.AdjustForCentreOfMass(1e-3);
        
        if (opts.PotFile != "")
            Energy.SetPotentialFromData(ref s, Table.Read(opts.PotFile));
        else
            Energy.CalculatePotentialShell(ref s, opts.Softening);
        
        // make sure energy is calculated for each particle
        double e_tot;
        Energy.CalculateTotalEnergy(ref s, out e_tot);
        
        Console.WriteLine(((EnergyParticle)s[0]).Energy);
        
        EnergyBin n = new EnergyBinDMDEDJ(new Bin.Function(EnergyBinDMDEDJ.E), new Bin.Function(J), null, 0, 0);
        n.NumBins = opts.NumBins;        
        n.Analyse(s);
        
        if (opts.Binary)
            BinaryOutput(n.Data, opts.OutFile);
        else if (opts.Maple)
            MapleOutput(n.Data, opts.OutFile);
        else if (opts.Gnuplot)
            GnuplotOutput(n.Data, opts.OutFile);
        else
            n.Data.Print(opts.OutFile);
    }
    
    public static double J(Particle p)
    {
        return Math.Sqrt(p.J2);
    }
    
    public static void GnuplotOutput(Table data, string file)
    {
        using (StreamWriter sw = new StreamWriter(File.Open(file, FileMode.Create)))
        {
            int n = (int)Math.Sqrt(data.Length);
            for (int i = 0; i < data.Length; i++)
            {
                if (i % n == 0 && i > 0)
                    sw.WriteLine();
                for (int j = 0; j < data.Columns; j++)
                {
                    double x = data[i, j];
                    sw.Write(String.Format("{0:f6} ", x));
                }
                if (i < data.Length - 1)
                    sw.WriteLine();
            }
        }
    }
    
    public static void MapleOutput(Table data, string file)
    {
        using (StreamWriter sw = new StreamWriter(File.Open(file, FileMode.Create)))
        {
            int n = (int)Math.Sqrt(data.Length);
            for (int i = 0; i < data.Length; i++)
            {
                if (i % n == 0 && i > 0)
                    sw.WriteLine();
                sw.Write(String.Format("{0} ", data.Z(i)));
            }
        }
    }
    
    public static void BinaryOutput(Table data, string file)
    {
        // no header, offset is zero
        // doubles because they're standard
        // use wip command: image filename`colsdrows
        using (BinaryWriter bw = new BinaryWriter(File.Open(file, FileMode.Create)))
        {
            int n = (int)Math.Sqrt(data.Length);
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    bw.Write(data[i + j * n, 2]);
            }
        }            
    }
}
